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Abstract 

Using straightforward linear algebra we derive response operators describing the impact 
of small perturbations to finite state Markov processes. The results can be used for studying 
empirically constructed - e.g. from observations or through coarse graining of model simula¬ 
tions - finite state approximation of statistical mechanical systems. Recent results concerning 
the convergence of the statistical properties of finite state Markov approximation of the full 
asymptotic dynamics on the SRB measure in the limit of finer and finer partitions of the phase 
space are suggestive of some degree of robustness of the obtained results in the case of Axiom 
A system. Our findings give closed formulas for the linear and nonlinear response theory at 
all orders of perturbation and provide matrix expressions that can be directly implemented in 
any coding language, plus providing bounds on the radius of convergence of the perturbative 
theory. In particular, we relate the convergence of the response theory to the rate of mixing of 
the unperturbed system. One can use the formulas obtained for finite state Markov processes 
to recover previous findings obtained on the response of continuous time Axiom A dynamical 
systems to perturbations, by considering the generator of time evolution for the measure and 
for the observables. A very basic, low-tech, and computationally cheap analysis of the re¬ 
sponse of the Lorenz ’63 model to perturbations provides rather encouraging results regarding 
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the possibility of using the approximate representation given by finite state Markov processes 
to compute the system’s response. 

1 Introduction 

1.1 A Brief Summary of Response Theory 

The development of methods for computing the response of a complex system to small pertur¬ 
bations affecting its dynamics is a subject of very active investigation in many fields of science 
and of technology. Statistical mechanics provide the tool for approaching such a problem 
through so-called response theories, which allow for evaluating the change in the properties of 
a system through suitably defined operators - response formulas - that factor in the statistical 
properties of the unperturbed system and the specihc nature of the perturbation one wants to 
study. 

One can see a response theory as a virtual experimental setting where one has at hand a 
given system, various measurement instruments, and a knob controlling the value of a param¬ 
eter, and knows how to relate the position of the knob with the reading of the instruments. 
In other terms, response theories provide the basis for understanding the outcome of exper¬ 
iments, and, not by chance, physical sciences have been at the forefront of the theoretical 
investigation in this direction. The monumental contribution by pQ has provided the basis and 
the explicit formulas needed for studying the impact of very general perturbations to statis¬ 
tical mechanical systems at equilibrium, as described by the canonical ensemble. The Kubo 
formulas are extremely useful for studying a large class of problems in e.g. transport, optics, 
and acoustics. A cornerstone of Kubo’s theory is the fluctuation-dissipation relation, which 
enables connecting - within linear approximation - the free fluctuations of the system to its 
response to perturbations. This property is closely related to the celebrated diffusion law for 
the brownian motion and has been recently extend to a fully nonlinear case [2]. Despite its 
obvious relevance, Kubo’s approach has been criticized for several reasons; 

• it is not physically consistent in treating the transition from equilibrium to non-equilibrium 
dynamics, because it studies the impact on equilibrium systems of perturbations that 
drive them near (but out of) equilibrium, but does not clarify how a new stationary state 
is reached and maintained, and at the same time it is not suited for studying the response 
to perturbations of non-equilibrium systems; 
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• it lacks mathematical rigour, as it is not clear which are the systems for which the 
response formulas apply, and why it should apply at all. 

In [3l U [5] it was clarified that it is possible to establish a rigorous response theory for Axiom A 
[6] continuous or discrete time dynamical systems. One obtains that the invariant SRB measure 
is smooth with respect to the parameter e that controls the strength of the perturbation 
changing the dynamics of the system from x = F{x.) to x = F(x) + eX(x), in the case of 
continuous time evolution, and from x^+i = F(xfc) to x^+i = F(xfc) + eX(xfc), in the discrete 
case. We continue our discussion in the continuous time case. 

We can introduce the unperturbed evolution operator Sq = exp(tF-), which moves forward 
in time any function of phase space 0(x) by an interval t according to the unperturbed 
dynamics, so that 0(x(t)) = S'qO(x( 0)), and its perturbed counterpart 5* = exp(f(F + eX)-), 
which instead describes the evolution in the perturbed system. 

We define /9o(dx) and /5e(dx) the invariant measures of the unperturbed and perturbed 
states, respectively. In particular, one obtains that the expectation value of sufficiently smooth 
observables 0(x) in the perturbed state can be expressed in the form: 

OO 

[0]. = |0|„ + 53eJi[0],, (1) 

i=i 

where [Q]^ = f i/f:(dx)Q(x) and [QJo = f iyo(dx)Q(x), while the various terms of the pertur¬ 
bative expansion can be written as: 

/ roo poo 

z.o(dx) / dti.../ dt^ASlF..Sl”-^ASl-Oix), (2) 

Jo Jo 

where A(») = X • V(»). In particular, the linear term can be written as: 

dtiA5*iO(x), (3) 

All terms S[0]j can be written as an expectation value on the unperturbed measure of a new 
observable expressed as a functional of the background vector field F, of the perturbative 
vector field X, and of the observable O. The somewhat surprising conclusion we draw is that 
the invariant measure of the system, despite being supported on a strange geometrical set, is 
differentiable with respect to e. Among the many merits of the Ruelle response theory, one 
can mention that a) it clarifies the mathematical framework needed for developing a response 
theory, whose main ingredient, roughly speaking, is the robustness deriving from having a 
uniformly hyperbolic dynamics on the attractor supporting an SRB measure; and b) it works 
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seamlessly, in principle, in equilibrium and non equilibrium statistical mechanical systems, 
reducing to Kubo’s formulas when considering the first scenario, if one assumes that statistical 
mechanical systems are Axiom A. Non-trivial implications of the nonequilibrium/equilibrium 
dichotomy regarding the validity of the fluctuation-dissipation relations are discussed in 130 

E], 

Of course, at this stage one needs to bridge the gap between mathematical formalism and 
physical meaningfulness. One manages to bring Ruellle’s formalism into the realm of applica¬ 
bility by adopting the chaotic hypothesis PE], which basically says that a high-dimensional 
chaotic physical system can be treated at all practical purposes as if it were Axiom A if we 
focus on macroscopic observables. The chaotic hypothesis is the generalisation of the ergodic 
hypothesis, and provides a firm background for translating the mathematical properties of 
Axiom A systems into physically meaningful statements. Clearly, the chaotic hypothesis ap¬ 
plies far from regimes of metastability and far from critical transitions, where entirely different 
phenomena appear. The chaotic hypothesis might also be practically problematic in the case 
one treats multiscale systems featuring many near-zero Lyapunov exponents; see discussion in 

m- 

Taking the point of view of the chaotic hypothesis, one has that, after transients have died 
out, nonequilibrium systems reach a nonequilibrium steady state (NESS) where the phase space 
is on the average contracting (with the rate of contraction corresponding, broadly speaking, 
to the entropy production of the system m), so that one can associate to the hyperbolic 
strange attractor supporting the invariant measure a Hausdorff dimension that is lower that 
the dimensionality of the phase space and, in general, not integer mm- 

The last piece of the puzzle one needs to lay in order to sort out the above-mentioned 
criticisms to Kubo’s theory relies on the physical interpretation of how a perturbed equilibrium 
system reaches a steady state. A convincing point of view on this relies on emphasizing the 
role of thermostats, which are large physical systems interacting with the system of interest 
in such a way to extract the excess of heat generated as result of the energy input due to the 
perturbation. Thermostats are also responsible for making it possible the set-up of stationarity 
in the case of forced and dissipative non equilibrium systems. An extensive treatment of the 
role of thermostats in equilibrium and nonequilibrium systems in the context of the chaotic 
hypothesis is given in m- We will not elaborate further on this aspect here. 
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1.2 Transfer Operator Approach 

One can point out that the formulas above describe the impact of and expressed in terms of 
expectation values of a generic observable O, whereas one might like to derive directly results 
for the impacts of the perturbations on the invariant measure. 

In O in [5] one constructs the response of the system to perturbations by following the 
changes in the individual trajectories and summing over the possible initial configurations 
distributed according to the unperturbed invariant measure. A different point of view on 
response theory focuses on studying the properties of the unperturbed and perturbed transfer 
operators and of their generators (see |14) for an introduction on these mathematical objects), 
through the construction of an appropriate framework of suitable (Banach) functional spaces 
where their actions are well defined, able to carefully treat the fundamental differences between 
the (smooth) unstable and (singular) stable manifolds of the Axiom A systems |T5l [161 El HB] • 

The evolution of the measure driven by the system x = F(x) up to time t > 0 starting from 
an initial condition at time t = 0 is described by the Perron-Frobenius transfer T* (see, e.g., 
[II]), so that p(x, t) = £*/?(x, 0). We have that the family of {C^}t>o forms a one-parameter 
semigroup, such that and £^ = 1. The Perron-Frobenius operator £* is the 

adjoint of the evolution operator 5* = (T*)"*", so that {S^O,p) = {0,C^p), where {f,g) is 
the action (computation of the expectation value) of the linear functional g (the probability 
measure) on the test function / (the observable). We have that = i^o > 0, meaning 
that the invariant measure is an eigenvector corresponding to unitary eigenvalue of the Perron- 
Frobenious operator. 

Assuming strong continuity and boundedness of the semigroup given by {£*}t>o, we can 
introduce the unperturbed Liouvillian operator L, which is the generator of the unperturbed 
Perron-Frobenius operator £* = exp(tL), and write the Liouville evolution equation for p(x, t) 
as follows HU: 

dtp = -V • (pF) = Lp (4) 

One immediately obtains that Lvq = 0. In general, the spectrum of L is complex and in a 
strip of finite width including and below the imaginary axis consists only of isolated eigen¬ 
values of finite multiplicity corresponding to the Ruelle-Pollicott resonances, while below such 
strip one finds the essential spectrum, which is responsible for the continuum of the power 
spectra of integrable observables. Furthermore, the presence of a unique SRB measure comes 
from the presence of a simple vanishing eigenvalue, while mixing properties result from the 
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absence of any other eigenvalue along the imaginary axis. The relevance of these properties 
for constructing a response theory are discussed in great detail in [mils]. In [20] it is argued, 
using mathematical considerations and examples of geophysical relevance, that the presence 
of Ruelle-Pollicott resonances having real part close to zero may lead to the presence of rough 
parameter dependence, as the smoothness of the response if lost. Additionally in [2T], it is 
shown, along similar lines, that the crisis of a very high-dimensional chaotic attractor near a 
critical transition - namely, of a climate model in the vicinity of the tipping point responsible 
for the transition between warm and snowball climate ISa EallSl ESI - can be detected and 
anticipated by looking at spectrum of the transfer operator. 

We then have that the presence of the e perturbation to the dynamics changes the Liouville 
equation as follows: 

dtp = -V • (pF) - eV • (pX) = L,p, (5) 

so that we can introduce the perturbed Perron-Frobenius operator C\ = exp(tLe), which 
pushes forward in time the measure according to the perturbed dynamics: p(x, t) = 0). 

Clearly, {S\0,p) = {0,C\p). Additionally, we have that = Pe Vt > 0 and = 0. 
While this approach is in some sense mathematically more problematic, because it is based 
on studying a partial differential equation instead of a finite dimensional dynamical system, it 
seems to provide a more comprehensive set of tools for studying the response of a system and 
relating it to its unperturbed fluctuations, see, e.p., [I5|, where Ruelle’s formulas are obtained 
along these lines. See also a comprehensive review given in where the applicability of the 
response theory beyond the case of Axiom A systems is discussed in detail.. 

One needs to emphasise that the transfer operator approach is more natural in all the 
cases when our interest focuses on studying the properties of the response of an ensemble 
of trajectories (initialised according to the unperturbed invariant measure) rather than on 
individual orbits of a system. 

Note that in some applications there is not an obvious separation between the two ap¬ 
proaches. Let’s take the problem of constructing climate projections through the use of (ex¬ 
tremely complex) numerical climate models, which is one of the core activities summarized 
in the IPCC reports [26]. Indeed, modelling centers are actively pursuing the preparation of 
multiple runs starting from an ensemble of initial conditions for a given scenario of forcing in 
order to estimate more accurately the uncertainties in the projections. Nonetheless, we will 
not experience an ensemble of realizations of the climatic evolution, but just one. 
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1.3 Computing the Response 

The analysis of high-dimensional complex system in terms of direct numerical simulation and of 
time series analysis suffers from the (almost) ubiquitous curse of dimensionality, which makes 
it hard to represent correctly the details of the dynamics because computational complexity 
explodes with the number of degrees of freedom. The construction of efficient and accurate 
algorithms for studying the response of a complex system to perturbations faces serious dif¬ 
ficulties. Let’s focus now on the linear case. Some previous studies has emphasised the need 
for treating separately the contributions to the response coming from short and long-time de¬ 
layed contributions in Eq. and have underlined the need for reducing the complexity of the 
invariant measure by adding in the background state some stochastic forcing, able to smooth 
our the singularity of the SRB measure [271 [2H] • 

A promising way to deal with the actual computation of the scalar product in Eq. [^is to use 
as time-dependent basis the covariant Lyapunov vectors IMIEU], which automatically separate 
the contributions to the response coming from the unstable, neutral, and stable directions. 
Tis clarifies that the convergence of the formula given in Eq. comes from the two distinct 
facts that a) perturbations along the stable directions naturally decay, and b) perturbations 
along the unstable directions grow in size, but are dominated by the loss of correlation dne to 
mixing. 

Recently, algorithms based upon adjoint methods have shown a good degree of accuracy 
and seem promising, even if scaling them up to high-dimensional systems has not been at¬ 
tempted yet (SUES]. A different approach to the problem has been proposed in IMllZlIMlESj, 
where, instead of trying to computing ab initio and directly the response given in Eq. the 
authors construct it a posteriori, probing the system with some test forcings and using the 
formal properties of the theory to be able to predict the response for new patterns of forcings. 
One can say that by studying the differential response to similar yet differently modulated 
perturbations, it is possible to derive the overall response properties of the system. 

1.4 This Paper 

Any numerical representation of a continuum system builds upon the need of discretizing the 
phase space and, in the case of time-continuous system, of time. 

In this case, we partition the phase space of the system in say n states 4>i,.. .(pn- In many 
cases, the states are constructed by discretizing the phase space in a grid of boxes, which 
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provide a (Galerkin) basis of orthogonal functions. We then construct an initial ensemble as 
defined by the occupancy Uq, ..., Uq of each of the 0i’s, i = 1,..., iV, so that 


«o = 


/ 


dxp(x,0)l((/>j) 


where l(^) is the characteristic function in the set A, and we want to approximate the evolution 


of such occupancies change with time, considering discrete time steps At, so that, to a good 
approximation, 



Moreover, in such a discrete representation, we have that the value of an observable O in the 
state (j)i is given by its average 


n _ J dxp{x,kAt)l{(j)i)0{x) 
Jdxp(x,feAt)l(</.,) 


( 6 ) 


Let’s emphasize that when analyzing virtually any sort of complex system, almost invariably 


one proposes a natural spatial and temporal cut-off, so that one on not in fact interested in 
really being able to compute the response of any possible observable defined at any possible 
spatial and temporal resolution, whereas meso- or macroscopic properties are relevant. Going 


again to the useful example of climate science, it is commonly regarded as a good and useful 
question to learn about the change in the surface temperature in response to climate forcing on 


a spatial scale corresponding to say a continent or a fraction thereof, and on a temporal scale 
of say one year. Nobody would find useful nor intelligent to study the surface temperature 
response over extremely small temporal and spatial scales. 

Empirically, using long numerical integrations and defining the set of hnite states (pi, i = 
1,..., A, we can construct the stochastic matrix Mij describing the probability of performing 
a transition from state (pi to state pj in a period of time At. The same operation can in 
principle be performed using experimental and observational data. A fundamental issue at the 
core of such procedure is whether for some dynamical systems in the limit of finer and finer 
partitions covering the phase space (actually, the attractor of the system) with A —)• oo one 
reconstructs the actual invariant measure of the original system. See in [36] a comprehensive 
discussion of such an issue, the so-called Ulam conjecture, and in m some extremely promising 
applications of finite state Markov processes for studying severely reduced representations of 
complex systems. 


Following the idea that the performing the discretization of the phase space amounts to 
adding a stochastic perturbation of the original dynamical systems, with intensity going to 



zero with the scale of the actual partitions, and exploiting the fact that the SRB measure can 
be constructed as zero-noise limit (with measure that is absolutely continuous with respect to 
Lebesgue) of the physical measure, in [38l |39] it has been proposed that the Ulam conjecture 
applies in the case of Axiom A systems, which are endowed with an SRB measure. The 
convergence in the case of Anosov diffeormorphism has indeed been proved provided one adds 
some noise of asymptotically vanishing intensity (through stronger than the noise induced 
by the partition itself) to the underlying dynamics |40] . Somehow this is not so surprising 
because by adding noise one introduces a cutoff below which partitions do indeed work. At 
any practical level, these results suggest that in the case of Axiom A system constructing 
finite state Markov processes using Ulam partitions can do a pretty good job in simulating the 
true dynamics, if one consider reasonably well-behaved, smooth observables as test functions. 
Nonetheless, one has to note that different choices for the partitions can lead to very different 
rates of convergence [36j . See also the discussion and the numerical example presented in |41j . 

Apart from the Ulam method, one can follow a mathematically more elegant but practically 
much harder way to construct finer and finer partitions. As well known, Axiom A systems 
possess Markov partitions, i.e. well-defined, metric independent, finite resolution represen¬ 
tations of the phase space that refine themselves with the dynamics mm- Such Markov 
partitions can be used to construct in the limit the actual SRB measure of the system, and, 
additionally, following [32], they provide a natural way to build finite Markov chains whose 
properties converge in the limit to those of the Perron-Frobenius operator of the system. 

Having a response formulas in the finite case has direct relevance for finite Markov chains 
and for interpreting the results of reduced models. Another good reason to construct a re¬ 
sponse theory in a finite state space has to do with the fact that the response operators for 
Axiom A systems introduced by Ruelle can be written as expectation value of certain observ¬ 
ables on the unperturbed SRB measure. Therefore, given what said above, one can hope to 
have convergence of the finite state reconstructed response operators to the corresponding true 
response operator in the limit of infinitely fine partitions of the dynamics. Actually, providing 
explicit formulas for the response operator for a finite state partition of a system the response 
operator and taking the limit for (suitably defined) finer and finer partitions could be inter¬ 
preted as a rigorous way for constructing the actual response on the asymptotic SRB measure. 
One needs to note - see discussion in Sects. 2T and[^- that special attention has to be paid 
when studying the convergence of such operators. 

In what follows, we present the derivation of the response formulas at all orders of pertur- 
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bations (as well as the full nonlinear versions) for finite state spaces of arbitrary size N. All 
expressions are given in terms of the transitions matrix of the unperturbed system, to its cor¬ 
rections to the perturbation, and of the parameter controlling the strength of the perturbation. 
The interest we see in the calculations we present below is mostly three-fold: 

• our results are obtained using basic linear algebra operations in finite dimensional spaces, 
which can used to interpret more complex operators acting on infinite dimensional spaces. 
It is also possible to use the finite dimensional expressions to derive, e.g., the the actual 
response operators for continuous time Axiom A dynamical systems; 

• we are able to derive an explicit expression for the a lower bound to for the radius of 
convergence of the perturbative theory, and relate it with the mixing properties of the 
unperturbed system. We also find a (very tentative) expression for such a lower bound 
in the case of continuous time case Axiom A dynamical systems; 

• our formulas can be translated into one-line commands in now widely available software 
tools like R, Octave, or MATLAB®. This might greatly facilitate the actual implementa¬ 
tion of response operators. In particular, we can say that our results provide a direct 
translation of the response theory into a readily implementable algorithms. 

The paper is organised as follows. In Sect. we introduce some notation and provide basic 
properties of ergodic finite state Markov chains, which can be taken as mathematical model 
on its own or as finite precision representation of ergodic (in this case. Axiom A) systems. 
We also show how it is possible to find an exact expression for the impact of a perturbation 
on the invariant measure of the Markov process and we study the radius of convergence of 
the perturbative expansion. In Sec. we rephrase our results in terms of observables, by 
constructing straightforward adjoint operators in finite dimensions. In Sect, [^we show how 
our findings agree with the response theory for continuous time systems when we suitably 
translate the matrix operations into operators. In Sect, [^we present a simple yet instructive 
investigation of the response of the Lorenz ’63 system |43j to perturbations using Ulam-like 
partitions and the formalism developed here. In Sect, [^we recapitulate and discuss our results. 
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2 Response Operators for Finite-state Markov Pro¬ 


cesses 

Let’s consider an ergodic Markov process with a finite number of states defined by the N- 
component vector u. We consider the inhnite Markov chain generated as uq, A^uq,. • • Af”uo, 
... where uq is the initial ensemble of states, and Aiij G is the stochastic transition 

matrix determining the probability of reaching the state i at step n if at step re — 1 we are in 

the state j. The process is taken to be stationary, so that M. does not change with re. We 
remind that A4 is such that 1 M-ij > 0 Vi, j = 1,..., V. 

The invariant measure is obtained by solving the eigenvalue problem 

Alu = Au, (7) 

and selecting the unique solution with eigenvalue A = 1. The corresponding (column) eigen¬ 
vector uij^ is the invariant measure of the system.We also remind that 

lim M"'z = aiui, Vz. (8) 

71^00 

where {Aj, uj} j = 1,..., re are the pairs of eigenvalues and eigenvectors of A4, where Ai = 1, 

|Aj| < 1 if j > 1, and z can be expressed as z = 

Our goal is to find a formula for expressing the change in the invariant measure resulting 
from perturbing the transition matrix A4 ^ Ai + em. 

We note that in order to preserve the Markov property of the system, rre obeys the fol¬ 
lowing constraint: ~ O’ = 0- Moreover, an additional 

constraint on em comes from the fact that all elements of At -|- em have to be positive. We 


define 

e+ =max|Vz,j G {1,..., 

e 

N},Aiij + emjj > 0, 

(9) 

and 

e_ = min Vz, j G {1,..., 

e 

N},Aiij -b emij > 0; 

(10) 

clearly, 

to have 

e_ < 0 < e+, and the perturbed matrix is a stochastic matrix Ve G [e_, 

some room for studying the impacts of perturbations, we require that 

e+]. In order 

e+ — e_ > 0. 


^Most commonly Markov chains are constructed using row vectors; we use column vectors because we find it easier 
to perform formal matrix manipulations and because we are closer to the formulation most commonly implemented 
in scientific software. 
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Such conditions show that, for a given M, it makes sense to consider only a specific class of 
perturbation matrices m. Let’s provide an example of an ill-chosen m: if Ai has two zero 
entries = 0 < 0, then we have e_ = 0 = e+. 

The new invariant measure is the unique solution to the eigenvalue problem: 

{M + em)u = Au, (11) 


with unitary eigenvalue. We dehne vi as the invariant measure of the perturbed system. 
Our goal is to express it as a function of Ad, m, e and u. This amounts to constructing a 
response theory. We first present the results of the explicit calculation, and then discuss issues 


of well-posedness of the problem and convergence of the procedure in Sect. 2.1 Let’s express 
vi = ui + e-^ wJ, so that we obtain: 


CO CO 

(Ad + em)(ui -|- e-^wj) = ui + e-^w^, (12) 

i=i 1=1 

Note that the first eigenvalue is not changed by the perturbation Ad —)• Ad -|- em, because also 
Ad + em is a stochastic matrix. Using the definition of ui we obtain a system of concatenated 
equations 


(1 — Ad)w^ = mui 
(1 — Ad)w" = 


Vn G N, n > 1 


We obtain 


= ikiui = (1 — Ad) ^mui 


n-l 


w" = TiW 


Given the recursive structure, we immediately derive the general formula: 


w" = = n (^(1 - -Ad) ^ mj ui. 

j=i 


where Concluding, we have that: 


oo n 


(13) 

(14) 

(15) 

(16) 

(17) 


Vi = Ui + ^ e"wn = Ui + ^ = ui + e" ((1 “ -Ad) ^ mj ui (18) 

n=l n=l n=l j=l 

which provides the formula we have been looking for. We note that the term responsible for 
the order of perturbation to the measure can be expressed as 


1 d^ 

hm — -—vi. 
e->o n! de"’ 


( 19 ) 
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Using the matrix identity (1 —M) ^ ~ ^ (1 “ -^) ^ m, we can also 

formally express the previous result as: 

vi = (1 — e'I'i)“^ui = (1 — e(l — (20) 

Using again the matrix identity (1 — M)~^ = previous expression can be 

rewritten as: 

OO 

vi = (1 — = (1 — (21) 

or 

OO n / OO \ 

vi=^i+E^"n E M^m ] ui (22) 

n=l j=l Vfc=0 / 

2.1 Well-posedness and Convergence 

In the previous equations, we have used somewhat carelessly the expression (1 — M)~^. Un¬ 
fortunately, the matrix 1 — A1 is not invertible, because all of its columns sum up to zero, 
or, alternatively, because we know that 1 is an eigenvalue of Ai. Nonetheless, the expression 
makes sense if we apply it to a vector belonging to span{u 2 ,..., Un}. We now want to prove 
that: 

Lemma 1 If Ai is a Markov transition matrix —)• with eigenvectors (ui,U 2 ,... ,Un), 

and corresponding eigenvalues (Ai = 1, A 2 , • • •, A,i), 1 > IA 2 I > ...|A„|, and m is a matrix 
matrix —)• M” such that = 0? then mz G span{u 2 ,..., Un} Vz G M"'. 

Proof Let’s consider the vector y = mz. Its component can be written as t/i = 
X]j=i Since mij = 0, we have that Y^=i = Yf=i Yf=i = 0- 

Let’s now consider the eigenvector u^. of Ai. We have Y]=i Since 

Y7=i-^i,j ~ taking the sum over the i components of the previous expression, we obtain: 
Yl=i Y!j=i -M-ijUk-j = Y%i '^k-,j = Afc Therefore, either \k = 1, or X)”=i Uk-,j = 0. 

We have that if /c > 1, Y^=i '^kj = 0. 

We conclude that y = mz G span{u 2 ,..., Un} Vz G 

Remark One needs note that finite numerical precision might cause troubles, so that one 
should be careful in eliminating any component along ui at each before applying X^^i AI-^. 
Note that we must use Xl^i expression for (1 — Ai)~^ in any code, because otherwise any 
software would give us automatically a NaN as error message. 
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Remark We wish to underline another method for avoiding the NaN problem mentioned 
above. Following |44j . we introduce the fundamental matrix of the Markov chain as 2^ = 
(1 — A4 + where is the limit matrix whose columns are all equal to u. One can 

show that 2 exists as the operation of inverse is well defined given the spectral properties of 
A4 — A4°° [38j- One can show that Ai°°mz = 0 Vz G Therfore, in all the previous Eqs. 

D 


16H22 


we can substitute (1 — M.) = Zm = YT=o{-^ ~ 


Let’s consider the problem of convergence of the expression in Eq. 18 We want to make 
sure that the norm of does not diverge, and use this to find a bound for the 

value of e. A simple way to approach this problem is to study the ratio of the norm of two 
consecutive terms in the previous series. Using Eqs. TSpU we have: 


\Wn\\l 


e"- i||u;n_i||i 


= e- 


l-M) ^ ^11(1 - A4) 


\Wn-l\\l 


\Wn-l\\l 


< e||(l - M) 


(23) 

(24) 

(25) 


where we use the submultiplicative property of the norm and we introduce a modified definition 
of the norm taking into account that the vector mv G span{u 2 ,..., Un}Vv G 


IIQIIt = 


sup 


IIQv||] 


r)espan{u 2 ,...,Un},||v||i=l I l''^N 1 


Using expression 25, we have that the perturbative expression converges if 


|e|(l - ||M||D"'IIHIi < 1 ^ kl < (26) 

The previous expression provides an explicit bound for our calculations. We note that €max is 
finite because of the restriction imposed in the definition of the norm || • ||*. Such a bound 
ensures also the invertibility of (1 — eTi)”^. From the previous result, we find the following 
bound for the first order correction to 


k'w^illi < 


e \m\\i 


i-\\M\\r 

so that ||m||i/(l — ||A4||*) can be though as a bound to the first order sensitivity of the measure 
to perturbations. 


Using expression 24 we can derive a more generous bound for e: 

1 


e ^ 


|m||i||(l — My 


> Cr; 


(27) 
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while ||m||i||(l —provides an additional (stricter) bound to the first order sensitivity. 
Note that in all the previous expressions we can substitute ||(1 — with ||2^||i. 

At this point, we wish to refer to previous results (see, e.g., [l5]) providing bounds for the 
norm of the difference between the perturbed and unperturbed invariant measure: 

,, ,, e||™,||i 

||wi - Willi < - JJT (28) 

where is the so-called ergodicity coefficient [l6] defined as: 

= ^sup||Af(ei -ej)||i 

with ej indicating the unit vector having 1 at the entry and zero elsewhere. We remind that 
ri(Af) is larger than any subdominant eigenvalue of A4, and 1/(1 — can be taken as a 

definition of conditioning number of Ai m- Clearly if T^i^) is close to 1, the bound given in 
Eq. [28| diverges. Note that 1/(1—r^(l)) is the bound to non-perturbative sensitivity mirroring 
the bound to the perturbative, linearized sensitivity given previously as 1/(1 — ||Af||)^). See 
also additional results presented in 


The sensitivity of the unperturbed measure to perturbations given in Eq. 28 can also be 


cast in terms the smallest possible value for constant controlling the rate of convergence 
of iterates Adoj, Ai^ei, ..., to ui, so that Vn G N+, Vz G 1,... A we have that ||Af"'ei — 

uilli < Cp"^, C > 1 [inillT]. The sensitivity diverges as pM approaches 1, i.e. when the 
unperturbed matrix has slow properties of convergence. 

While the quantities ||Af||^, and pM are indeed different, they all point to the fact 

that if the mixing rate of the unperturbed matrix Ai is slow - so that such quantities are close 
to 1 (so that 11(1 — Af)“^||| and \\Z\\i are very large) - then the sensitivity of the measure 
to perturbations is high. See in [JOj a discussion of the link between slow mixing of a system 
and the presence of rough parameter dependence in its response to perturbations, with some 
examples of applications in a geophysical context. 

Bringing together the results presented in Eqs. [9 10 and in Eq. 27, we conclude that Eqs. 
[T8]|^ provide the exact expression for the invariant measure of the stochastic matrix Ai + em 


Ve G {[-e. 


max'! ^max 


]n[e_,e+]}. 


3 Response Theory for Observables 

Let’s now look at the problem in terms of impact of the perturbation m on the expectation value 
of observables. Observables live in the dual space of the densities, and, given our convention, 
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they are row vectors. They are approximated as having a constant value within each cell of the 
chosen partition of the phase space. The expectation value of the observable vr with respect to 
a measure w can be written as (vr, w), where (•,•) denotes the scalar product. By definition, 
we have that (vr, ^w) = (A''~vr, w), where indicates the transpose (and adjoint, because we 
are studying real functions) of A. 

Let’s look at the change in the expectation value of the observable vr as a result of At —>■ 
Ai + em. We can write: 


OO 

(vr, vi) = [vr], =[vr]o + ^ e"'5[vr]„ (29) 

n=l 

OO 

= (vr,ui) + ^e'^(T([vr,Ui) (30) 

71=1 

OO 

= (vr,ui)+ ^€^^((^^7) vr,ui), (31) 

n=l 


where [vrjo = (vr, ui) is the expectation value of vr in the unperturbed system, [vr], = (vr, vi) is 
the expectation value of vr in the perturbed system, <5[vr]n is the order perturbation, which 
can be expressed as 

(5[vr]„ = liin^^(vr,vi). (32) 

e^o n\ 

Moreover, is the order adjoint response operator, acting on the observables, which can 
be written as: 

n / OO \ 

j=l \k=l / 


We can also wrote Eq. 31 


as: 


e^7) ^7r,ui) 

(34) 

/ OO \ T 


em^ )-V,ui) 

(35) 

\fe=i / 


em'''(l — Al''')“^)“^vr, ui). 

(36) 


where the last two expressions provide the nonperturbative formulas. 


Remark Equations and 31 provide at all orders the response formulas for the discrete 


Markov process studied here. If we are constructing empirically the discrete phase space, we 
expect that different choices of the partitions, corresponding to different approximate repre¬ 
sentations of the full dynamics, will deliver different results in terms of response. Hence, our 
results can be model dependent, which is acceptable, as we are starting from a subjective choice 
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on the way we approximate the phase space. In fact, one can empirically test the robustness of 
the obtained results against a set of given criteria by comparing whether the perturbations to 
a certain set of relevant observables weakly depend on the specific partition used. We present 
a very preliminary (and encouraging) numerical study performed on the Lorenz ’63 model |43] 
later in Sect. [H 

Moreover, as discussed in Sect. |1.4[ if we construct finer and finer partitions of for studying 
the response of systems whose unperturbed dynamics features an SRB invariant measure (most 
notably in the case of Axiom A systems), and indeed if we follow the self-refining Markov 
partitions of the dynamics, our results should converge to the exact response theory built 
upon the true SRB measure. 


One needs to note that Eq. [27| gives an estimate of the largest possible value of e for a given 
partition, but we are are not sure whether the minimum over all the finer and finer partitions 
of is positive - this corresponds to imposing the uniform - in - bound on the norm of 


11(1 - My -|q or ||.i-||i. 


or ||Z| 

In [3H] it is shown that convergence of the finite state measure constructed using the 
Ulam method to the actual SRB measure is realized when ||2^||i grows asymptotically not 
faster than log AI, where N is the number of states. The requirement we seem to have here for 
applying response theory here is unavoidably stricter because computing the response entails 
considering the expectation value of not necessarily well behaved observables, constructed 
through nontrivial operations of differentiation of the actual observables of which we want 
to study the sensitivity to perturbations, see Eq. and laiiiis]. This essential difficulty is 
exactly what motivates the point of view discussed in [13119], where a delicate analysis of the 
relationship between tangent space of the unperturbed dynamics, the perturbation flow, and 
of the observable allow to set up a robust framework for the response theory. 

Similarly, in our case, making the response theory work at practical level means hav¬ 
ing/choosing m and u in such a way that ||(1 — A4)“^||^ or ||2^||i grossly overestimates in 


terms of norm the effect of applying {1 — M) ^ or equivalently Z in, e.g., Eq, 


22 


Addi¬ 


tionally, a suitable choice of the observable vr can help avoiding potential singularities in Eq. 
^6] In other terms, response theory can work much more easily once we get rid of or cure 
pathological cases. 
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4 Towards Continuous Time Dynamical Systems 


We want to rephrase the previous results in the context of continuous time dynamical systems 
and derive some formulas previously presented in the literature concerning Axiom A systems. 
We coonsider a time continuous dynamical system of the form x = F (x) and study its response 
to the perturbation F(x) —)■ F(x) + eX(x). Correspondingly, as a result of the perturbation, 
the original invariant measure z^(dx) is changed into /i(dx). The Liouville equation describing 
the evolution of a given initial density of states p(x) for the unperturbed system can be written 
as 

dtp{x, t) = -V ■ (F(x)/9(x, t)) ■ (37) 

considering two instants of time separated by a small time interval dt, we have: 

/?(x, t + dt) = /9(x, t) — dtV • (F(x)/9(x, t)) = Mp(x, t) 

M = l + dtJ^ F=-V ■ (F(x).). (38) 

We understand that Ai is in this context the unperturbed Perron-Frobenius operator 
pushing forward the measure p from t to t + dt. When looking at the perturbed flow we have: 

p{x, t + dt) = p{x, t) — dtV • (F(x)p(x, t)) — dfeV • (X(x)/?(x, t)) 

= {M + em)p{x,t), ■ (39) 


where 


m = dtX A = -V • (X(x)*) 


(40) 


In this case, starting from Eq. and considering that no normalization is applied to the 
perturbation operator, it is possible to propose a definition of continuous time 


dynamics taking inspiration from Eq. 27 


1 


TT—1 II*’ 

Wb 


||A||b 

such that the perturbative expansion converges if e < e 


(41) 


,, where || • ||g describes the norm 


of the operator in the appropriate Banach space B it belongs to, while || • ||g is such that 
the computation of the norm excludes the SRB measure. Note that e^ax is finite if both 
||A||g|| and are finite. This expression is admittedly tentative. As mentioned before, 

the problem of selecting appropriate functional spaces for constructing the response theory 
for Axiom A systems along the lines of studying the perturbations to the transfer operator 
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requires a careful construction of suitable Banach spaces and of the related metrics [El [HI [Is] 
and is beyond the scope of this paper0 


4.1 Linear response 


We now want to derive the Ruelle response formulas for computing the linear correction to 
the invariant measure resulting from the perturbation. We write 

OO 

z^e(dx) = 1 ^ 0 (dx) + ^ e"'r'n(dx), (42) 

n=l 


where n indicates the order of perturbation. Let’s first go back to the first order term in Eq. 


M 


ew^ = e'i'iui = I 


mui. 


(43) 


^fc=i 


Each term of the form pushes forward up to time tk = kxdt what is positioned to its right. 
Summing over k in, in fact, amounts to looking forward in time. If we insert the definition of 
m given above, we get the integrating factor df, so that we obtain the following expression: 


j-OO 

z^i(dx) = — dtV ■ (X(x(f))i/(dx)), (44) 

Jo 

where the evolution takes place according to the unperturbed system, and we have used the 
invariance of j/(dx) with respect to such an evolution law. 

By going into the dual space of the observables, we have that the change in the value of 
an observable 0(x) from time t to time t + dt in the unperturbed system can be written as: 

^0(x(f)) = F(x) • VO(x(t)), (45) 

so that 

0(x(f + dt)) = 0(x(f)) + dfF(x) • VO(x(f)) = Al'''0(x(f)). (46) 

where the operator A4~^ = 1 + dtT^ = 1 + dfF(x) • V(»). Along the same lines, one derives 
that the perturbation operator mJ acting on the observable can be written as mJ = dtX~^ = 
dtX(x) • V(»). Furthermore, we introduce the following expansion for the expectation value 
of 0(x): 

OO 

n=\ 

^Following m, one might tentatively consider the norms of the operator acting between the Banach spaces B 2 ^q 
and Bi, 5 +i. 
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where [0]e is the expectation value in the perturbed system, [O]o is the unperturbed expecta¬ 
tion value, and the corrections are included in the summation. 

Applying this expression to the first order term in Eq. 3T]|33 


6[Tr]i = = e( j m’’' I j j 7r,ui). 


(48) 


\k=l 


we get: 


/ poo p poo 

i^(dx) J dtX(x) • VO(x(t)) = e J z^(dx) J dtA5QO(x) (49) 

which is exactly the original version of Ruelle’s linear response formula given in Eq. 

One needs to note that what in Ruelle’s formulation is causality (time integration in the 
response starts from 0), in the context of the Markov matrices formalism followed here comes 
from the algebraic expansion of (1—A4)“^. The issues of convergence mentioned in the original 
paper by Ruelle can be translated in the rate of mixing of the system as determined by the 


properties of A4 discussed in Sect. 2.1 


4.2 Higher order terms 

We can repeat the same construction to derive the higher order perturbation terms in the 
case of the continuous time dynamical systems. Inserting in Eqs. T5]|T6 the expression 38 for 


Ai and expression 39 for m, we obtain for the second order the following expression for the 
perturbation to the invariant density: 


i/ 2 (dx) = / dti / dt 2 V • (X(x(ti)Vxpi) • (X(x(fi + ^ 2 ))) i^(dx), 


(50) 


/o JO 

while the expression for the order correction reads like 

poo poo 

Vn{dx) = (-l)’^e’^ / dti... / dtnV ■ (X(x(ti)... ' (X(x(ti + ... tn))) z^(dx), 

Jo Jo 

(51) 

Considering the adjoint problem and computing the higher order corrections to the expectation 
value of the observable O, we derive the general response formula proposed by Ruelle 


6[0]n = e^ i^(dx) / dti... / dtnA5*i ... 5*-iA5*"0(x), 

J Jo Jo 


(52) 


as 


reported in Eq. [2j 
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5 A very basic numerical experiment 


In order to make a (very) preliminary assessment of the potential of some of the ideas presented 
in this paper, we have focused on investigating some properties of the celebrated Lorenz ’63 
system [33]: 

X = a(y — x) 

y= x{p-z)-y (53) 

i = xy — (3z 


where we have chosen the standard value for the parameters a = 10, p = 28, and (3 = 8/3. 
We remark that such a system is not an Axiom A, but instead a singular hyperbolic system 
m, which possesses a chaotic attractor and an invariant SRB measure [52]. In a previous 
publication [33], we have performed an analysis of the linear and nonlinear response of the 
Lorenz ’63 to perturbations, extending a previous investigation by Reick [53] . which makes 
us confident that response theory can be safely applied at all practical purposes also in this 
case. We consider the special case of time-indepedent perturbations to the dynamics resulting 
from substituting p —)• p + e in Eq. 53, so that the perturbation flow can be written as 
eX(x) = [0 ex 0]^. 

We have then identified a 3-dimensional box B containing the attractor, defined as B = 
{{x,y,z) G TZ^\x G [—20,20], y G [—30,30], z G [—0,50]}, and subdivided it, a la Ulam, in 
smaller boxes of size using a regularly spaced cartesian grid. We have considered partitions 
obtained using small boxes with linear dimension given by dx = 2 x j, dy = 3 x j, and 
dz = 2.5 X j, along the three directions, with j = 1,2,4, see Fig. [Tj This amounts to 
partitioning B into 8000/smaller boxes. Note that our construction delivers a much lower 
resolution with respect to what used in, e.g., [53] . 

We run the model with standard values of the parameters choosing as initial condition 
[1 1 l]"*^ (in fact, given the global attractivity and ergodicity of the Lorenz attractor, any 

initial condition can be chosen), and, after discarding a transient of 1000 time units, which 
brings us safely into the asymptotic regime, we run the model for 50000 time units with a simple 
Runge-Kutta 4*^ order adaptive scheme and obtain the output with time step of 0.001 time 
units. This takes less than 10 minutes in a today’s commercial laptop with standard specifics 
using MATLAB®. We present results at such a low level of sophistication in order to clarify that 
the appracch proposed here is rather robust and of relatively simple implementation. 
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Figure 1: Attractor of the Lorenz ’63 system with indication of the cartesian grids used for con¬ 
structing the partitions of its phase space. See text. 

As the box-counting dimension or capacity of the attractor of the model given in Eq. 
is do ~ 2.05, we expect that the number of boxes k = 1,...,A^ needed to cover the 
attractor decreases oc We obtain a slightly lower exponent ~ 1.9, which is perfectly 

acceptable as we are far from the asymptotic regime where the scaling given by do is realized. 

For each value of j, the boxes B^. dehne the discrete states k = 1,..., N^. By counting 
the number of times the trajectory is included in each state and normalizing we derive 
experimentally the asymptotic normalized occupancies Instead, by tracking the transitions 
between the various discrete states, we construct the estimate of the stochastic transition 
matrix Afp,g describing the probability that the state (j)q makes a transition to the state cjyp in 
one time step. By hnding the eigenvector corresponding to the unique unitary eigenvalue of 
■Mp,q: we hnd the invariant measure, which agrees up to very high precision with the empirical 
occupancy rate Uk computed from the trajectory. As a hrst step, we evaluate the expectation 
values of four meaningful observables given by and z, as obtained from the time 

integration of the Lorenz model and from its discrete representation in terms of Markov chain. 

Table shows that the agreement is rather good even when extremely coarse resolution is 


22 






Table 1: Expectation value of the observables z"^, and z and their linear response with respect 

to the perturbation p —?• p + e. The hrst row refers to the integration of the Lorenz model given 


in Eq, The other rows refer to the empirical discrete Markov chain constructed using boxes of 
different sizes. refers to the number of states. The linear response of the observables dehned in 


Eq. 32 has been obtained using Eq. 48 The derivative with respect to e is estimated using hnite 
differences with e = 0.1. See text. 



{x^) 

(y^) 


(z) 


s[y% 


6[z]i 

Lorenz ’63 Model 

62.9 

81.2 

630.0 

25.6 

2.8 

3.7 

50.3 

1.01 

MC, 3 = 1, = 770 

63.2 

82.0 

630.5 

23.6 

2.9 

3.8 

50.3 

1.01 

MC, 3 = 2,Nl = 205 

64.3 

84.2 

632.2 

23.6 

3.0 

3.5 

49.7 

1.02 

MC, 3 = 4, = 56 

71.3 

84.8 

637.5 

23.5 

2.9 

3.9 

50.1 

1.02 


used. 

We then show how to compute the response of the system to the perturbation due to the 
introduction of the vector field eX(x). We keep in mind that when continuous time dynamics 
is considered, there is a very simple linear relation between the perturbation flow and the 


corresponding perturbation to the Perron-Frobenius operator, see Eqs. ^ 40 


Therefore, we repeat the the steps described above for the e—perturbed flow (we choose 
e = 0.1 in order to be on the safe side in terms of convergence), compute the new stochastic 
transition matrices -Mp^q, and derive the perturbation matrices errip^q = Aip% — ■Mp,q- Once 
mp^q and Aip,q are known, we can use them to compute the response of the systems at all 


orders of nonlinearity using Eqs. 22 and 36 One needs to note that because of the non¬ 
infinite integration time considered, of the non-infinitesimal perturbation applied, and of the 
somewhat arbitrary choice of the boxes, it can happen that the original and perturbed flow may 
be characterized by a different number of discrete states. We have observed such a difference 
only in the case j = 1, involving one single extra state for the perturbed flow, with normalized 
relative occupancy (< 10“®). This problem can be easily sorted out by imposing a cutoff and 
removing from the the discrete description all states with very low. 

As discussed above, one needs to test accurately the well-posedness and convergence of 
the expansion in order to be sure to obtain meaningful results. This is not our goal at this 
stage for such a preliminary numerical test of our results. Therefore, we limit ourselves to 


23 


























the less ambitious yet interesting goal of computing the linear response defined in Eq. 32 for 
the observables indicated above, using Eq. 48 The results are reported in Table and seem 
very encouraging. We have that the results are very stable with respect to changes in the 
resolution of the boxes, and agree to a high degree of precision with the results one obtains by 
empirically evaluating the sensitivity of the observables with respect to the introduction of the 
perturbation flow using two integrations, as well, in the case of the z observable, with what 
reported in |33j . We note that the results are virtually unchanged if one uses instead of the high 
resolution time series with time step of 0.001 time units sparser observations corresponding 
to, e.g. a time step of 0.01 time units. Obviously, using a time resolution lower by a factor 
of s with respect to what considered here, one derives by tracking the transitions a stochastic 
transition matrix corresponding to the power of the one obtained at higher resolution. 
This does not affect the results as long as the sampling is much higher than the characteristic 
time scale of the system, which can be approximated in 1/Ai 1.1 time units, where Ai is 
the positive Lyapunov exponent of the system. On longer time scales, instead, the stochastic 
matrix is quasi-degenerate, with all columns almost equal to the invariant measure 


6 Conclusions 

Taking the point of view of finite state Markov systems, we have been able to construct a 
perturbation theory for studying the impact of small perturbations to the background dynam¬ 
ics. While previous approaches focus on the constructing a theory able to account for the 
effect of adding small perturbations to the baseline flow, we focus on computing the change 
in the invariant measure and for the change in the expectation values of general observables 
(one problem being the adjoint of the other) occurring when the Markov transition matrix 
M ^ M + em. 

The perturbation term em has to be such that all the columns of the new stochastic matrix 
sum up to 1 and all entries are positive. All of our findings are obtained with rather simple 
linear algebra manipulations and using basic properties of the stochastic matrices. We can 
express the response as a perturbation series or, after suitable resummation, using compact 
exact formulas. We are also able to assess the convergence properties of the response theory 
by defining a value such that if |e| < the perturbative expansion converges. We 
have that the stronger is the mixing of the unperturbed system, the larger is the value of emax- 
These findings match well with previous results providing upper bounds to the sensitivity of 
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stochastic matrices to perturbations. 

Our results provide a direct algorithmic method for studying the response to perturbations 
for finite state Markov processes and have the advantage of allowing for an immediate and 
practical change of point of view between response theory seen in terms of changes of the 
invariant measure or in terms of changes in the expectation values of observables, by simply 
computing the transpose of the resulting finite dimensional linear operators. Our findings give 
closed formulas for the linear and nonlinear response theory at all orders of perturbations 
through explicit matrix expressions that can be directly implemented in any coding language. 

We can use our formulas to study the response to perturbations of finite state Markov 
processes constructed in order to have a simplified and treatable picture of a complex sys¬ 
tem. Given two different state spaces constructed using different finite partitions covering the 
attractor of the system, we cannot expect to obtain the same results for the change in the 
expectation value of a given observables. The results might indeed be model dependent, but 
this is the obvious price one has to pay because of the subjective choice of the reduced state 
space. An assessment of the robustness of the obtained results is key to applying our methods 
in the context of reduced models. Nonetheless, the extremely unsophisticated numerical study 
reported here on the Lorenz’63 model is quite encouraging at this regard, even if test should 
be made on much higher dimensional models. 

If the underlying dynamics is Axiom A (or Axiom A equivalent, as in the cases where 
the chaotic hypothesis applies), one can impose conditions such that the response operators 
constructed using finer and finer partitions converge to to the actual corresponding response 
operators constructed on the SRB measure. The conditions are stricter than what needed 
in order to have convergence of the unperturbed measure, the basic reason being that Ruelle 
response operators correspond to nontrivial observables. 

Our results can be thought as intermediate steps at finite precision leading to the correct 
response formulas in the limit. One needs to add as a caveat that going from finite state 
to functional spaces is far from trivial and requires a high degree of mathematical precision, 
which is beyond the scopes of this paper. Nonetheless, the finite construction proposed here 
seems to somehow point at why some important mathematical issues emerge when the Perron- 
Frobenius operator formalism is considered {e.g. selection of suitable norms for vectors and 
linear operators, definition of specific functional spaces for the observables). 

Interestingly, we can use the formulas obtained for finite state Markov processes to study 
the impact of perturbations to continuous time dynamical systems, after making a suitable 
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identification between the considered transition matrices and the evolution operators for mea¬ 
sures and observables. This operation is straightforward because there is a simple linear exact 
relation between the perturbation in the vector flow of the dynamical system and the pertur¬ 
bation in the Perron-Frobenius operator when infinitesimal time intervals are considered. As 
a result, we are able to derive in a very simple way previous formulas obtained studying the 
perturbations to the transfer operator as well as the original expressions proposed by Ruelle 
for the linear and higher order perturbations in the expectation values of observables. Using 
the results obtained in the finite state case, we propose a formula for the radius of expansion 
of the perturbative theory. 

One can envision that in the case the underlying dynamics is discrete, there is not such 
a one-to-one correspondence between perturbations to the vector field and perturbations to 
the Markov transition matrix. This can be easily checked when constructing the perturbed 
Perron-Frobenius operator resulting from adding a e correction to the vector field, which results 
into changes in the Perron-Frobenius operator at all orders in e. Therefore, the perturbative 
expansion is different in the two cases. Agreement is instead found in the limit e —)• 0, or, 
more practically, when we retain only the linear terms in e perturbative expansion, i.e. when 
aiming only at the linear response function. 

Future investigations will try, on the one side, to have a sharper mathematical look at 
the problem of going from hnite to infinitely small partitions of the phase space, and, on the 
other side, to delve in the numerical study of the effectiveness and efficiency of the proposed 
tools. Apart from testing the results on specific hnite state Markov systems, we will test how 
robust the proposed methods are when studying hnite state Markov processes that have been 
empirically constructed from time series of observations or of numerical simulations of high¬ 
dimensional complex systems. One may be led to hoping that it could be possible to have 
an accurate representation of the response of a high dimensional system to perturbations by 
constructing a smart hnite state model well suited to studying specihc observables of interest. 
Of course, in order to deal with the curse of dimensionality, one would like to be able to 
go beyond the Ulam method and deal with hnite partition of reduced phase spaces where 
projection is applied on many or even most dimensions. 

Our formulas may address the now long-standing problem of constructing suitable algo¬ 
rithms for studying the response of chaotic systems to perturbations. It is extremely hard to 
construct an algorithm for computing the (linear) response theory directly on the how, because 
serious problems emerge when considering the contributions coming from the unstable direc- 
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tions in the tangent space. This might have great relevance for studying problems, like climate 
dynamics, where a direct construction of the response operator is especially challenging and 
slightly indirect methods have to be used [Mj and a lot of effort has been devoted to defining 
the so-called atmospheric regimes and predicting their response to forcings |55j . 
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